home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / RAN4.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  1KB  |  56 lines
  1. FUNCTION ran4(VAR idum: integer): real;
  2. (* Programs using routine RAN4 must define the global variables
  3. TYPE
  4.    gl64array = ARRAY [1..64] OF integer;
  5.    gl65reals = ARRAY [1..65] OF real;
  6. VAR
  7.    glnewkey: integer;
  8.    glinp,glkey: gl64array;
  9.    glpow: gl65reals;
  10. in the main routine. The initialization block, IF (idum < 0), has been
  11. written in real arithmetic to avoid overflow on machines with 2-byte
  12. integers. With 4-byte integers, this block can be simplified with MOD
  13. and DIV. *)
  14. CONST
  15.    im=11979;
  16.    rm=11979.0;
  17.    a=430.0;
  18.    c=2531.0;
  19.    nacc=24;
  20. VAR
  21.    isav,j: integer;
  22.    jot: gl64array;
  23.    r4,dum: real;
  24. BEGIN
  25.    IF (idum < 0) THEN BEGIN
  26.       dum := idum MOD im;
  27.       IF (dum < 0.0) THEN dum := dum+rm;
  28.       glpow[1] := 0.5;
  29.       FOR j := 1 TO 64 DO BEGIN
  30.          dum := dum*a+c;
  31.          dum := dum-rm*trunc(dum/rm);
  32.          glkey[j] := trunc(2.0*dum/rm);
  33.          glinp[j] := trunc(4.0*dum/rm) MOD 2;
  34.          glpow[j+1] := 0.5*glpow[j]
  35.       END;
  36.       idum := round(dum);
  37.       glnewkey := 1
  38.    END;
  39.    isav := glinp[64];
  40.    IF (isav <> 0) THEN BEGIN
  41.       glinp[4] := 1-glinp[4];
  42.       glinp[3] := 1-glinp[3];
  43.       glinp[1] := 1-glinp[1]
  44.    END;
  45.    FOR j := 64 DOWNTO 2 DO BEGIN
  46.       glinp[j] := glinp[j-1]
  47.    END;
  48.    glinp[1] := isav;
  49.    des(glinp,glkey,glnewkey,0,jot);
  50.    r4 := 0;
  51.    FOR j := 1 TO nacc DO BEGIN
  52.       IF (jot[j] <> 0) THEN r4 := r4+glpow[j]
  53.    END;
  54.    ran4 := r4
  55. END;
  56.